#!/usr/bin/perl -w
use strict;
use Bio::SeqIO;
use Statistics::Descriptive;
#Copyright 2011. Jason Weirather.
#This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License v3 or any later version.
#get the names of all the fastas
my $outputname = shift @ARGV; 
chomp(my @names = `ls $outputname.*.fasta`);
#my $base_name_for_hmm = $outputname;
#new temp for this
my $tran = rand();
my $tname = 1;
if($tran =~ /0\.(\S+)$/) {$tname =$1; }
$tname = "temp_" . $tname;
foreach my $mcl_cluster (@names) {
  #print "$mcl_cluster\n";
  #print "aligning $mcl_cluster   ...\n";
  open (OTF,">$tname.fasta"); #sterilize labels
  my $in  = Bio::SeqIO->new(-file => "$mcl_cluster" ,
                           -format => 'fasta');
  my %lengths;
  while ( my $seq = $in->next_seq() ) { 
    my $id = $seq->id;
    $lengths{$id} = length($seq->seq);
    if($id=~/^(\S+)\s+/) { $id = $1; }
    print OTF ">$id\n" . $seq->seq . "\n";
  }
  close OTF;
  # this is where we add some parameters to
  # make the models a little sharper by dropping sequences that are on the
  # extreme high or low side of length
  my $numstdev = 2; # number of standard devations beyond the median to allow lengths for models
  my $stat = Statistics::Descriptive::Full->new();
  $stat->add_data(values(%lengths));
  my $avg_len = $stat->median();
  my $stdev = $stat->standard_deviation();
  my $maxlen = $avg_len + $numstdev*$stdev;
  my $minlen = $avg_len - $numstdev*$stdev;
  # for these fastas we will prune them to eliminate absurd lengths more than 2 (or the set number) over or under the standard devation.
  # for now we will only limit it on the high side, not the low side.
  open (OTF,">$tname.fasta2"); #sterilize labels
  my $in  = Bio::SeqIO->new(-file => "$mcl_cluster" ,
                           -format => 'fasta');
  my %lengths;
  while ( my $seq = $in->next_seq() ) { 
    my $id = $seq->id;
    if($id=~/^(\S+)\s+/) { $id = $1; }
    if(length($seq->seq) <= $maxlen) { print OTF ">$id\n" . $seq->seq . "\n"; }
  }
  close OTF;

  

  `clustalw -infile=$tname.fasta2 -outfile=$tname.aln -quiet -align -quicktree`;
  unlink("$tname.fasta2");
  #print "finished alignments\n";
  `./clustalTOselex.pl $tname.aln $mcl_cluster.selex`;
  `./selexTOstalkholm.pl $mcl_cluster.selex $mcl_cluster.stalkholm`;
  #print "finished converting to selex\n";
  unlink("$tname.aln");
  unlink("$tname.dnd");
  unlink("$tname.fasta"); #unneeded copy used to make clustalw work
  my $filesize = -s "$mcl_cluster.selex";
  if($filesize > 0) { 
    my $hmmbasename = $mcl_cluster;
    my $buildcommand = "hmmbuild -n $hmmbasename $mcl_cluster.hmm $mcl_cluster.stalkholm";
    `$buildcommand`;
  } else { 
    unlink("$mcl_cluster.selex"); 
    unlink("$mcl_cluster.stalkholm");
  }
  #unlink("$mcl_cluster.selex.temp");
  #`cat $mcl_cluster.hmm.temp >> $tablename.hmm`;
}
`zip -r $outputname-fasta.zip $outputname*.fasta`;
`zip -r $outputname-selex.zip $outputname*.selex`;
`zip -r $outputname-hmm.zip $outputname*.hmm`;
foreach my $mcl_cluster (@names) {
  unlink("$mcl_cluster");
  unlink("$mcl_cluster.selex");
  unlink("$mcl_cluster.stalkholm");
  unlink("$mcl_cluster.hmm");
}
